Lab 6 - Spatial Data

Goals:

Read in spatial data (shapefiles)
Simplify polygons (st_simplify)
Add and transform projection (st_transform and st_crs)
Create several finalized-ish maps in R
Join spatial data (st_join)
Find spatial intersections (st_intersect)
Interactive plots with tmap
Updating color schemes, base layers, etc.
Plotting spatial lines
Creating sf data from lat/lon coordinates

Awesome resource: Geocomputation in R by Robin Lovelace, available online: https://geocompr.robinlovelace.net/

library(tidyverse)

## ── Attaching packages ──────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──

## ✔ ggplot2 3.1.0     ✔ purrr   0.3.0
## ✔ tibble  2.0.1     ✔ dplyr   0.7.8
## ✔ tidyr   0.8.2     ✔ stringr 1.3.1
## ✔ readr   1.3.1     ✔ forcats 0.3.0

## Warning: package 'tibble' was built under R version 3.5.2

## ── Conflicts ─────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()

library(sf)

## Linking to GEOS 3.6.1, GDAL 2.1.3, PROJ 4.9.3

library(tmap)

## Warning: package 'tmap' was built under R version 3.5.2

library(leaflet)
library(ggrepel)
library(ggspatial)
library(RColorBrewer)
library(raster)

## Warning: package 'raster' was built under R version 3.5.2

## Loading required package: sp

## 
## Attaching package: 'raster'

## The following object is masked from 'package:dplyr':
## 
##     select

## The following object is masked from 'package:tidyr':
## 
##     extract

#preventing conflict in commands between packages
select <- dplyr::select
filter <- dplyr::filter

Useful information on file types (from gisgeography.com):

.shp is a mandatory Esri file that gives features their geometry. Every shapefile has its own .shp file that represent spatial vector data. For example, it could be points, lines and polygons in a map.
.shx are mandatory Esri and AutoCAD shape index position. This type of file is used to search forward and backwards.
.dbf is a standard database file used to store attribute data and object IDs. A .dbf file is mandatory for shape files. You can open .DBF files in Microsoft Access or Excel.
.prj is an optional file that contains the metadata associated with the shapefiles coordinate and projection system. If this file does not exist, you will get the error “unknown coordinate system”. If you want to fix this error, you have to use the “define projection” tool which generates .prj files.
.xml file types contains the metadata associated with the shapefile. If you delete this file, you essentially delete your metadata. You can open and edit this optional file type (.xml) in any text editor.
.sbn is an optional spatial index file that optimizes spatial queries. This file type is saved together with a .sbx file. These two files make up a shape index to speed up spatial queries.
.sbx are similar to .sbn files in which they speed up loading times. It works with .sbn files to optimize spatial queries. We tested .sbn and .sbx extensions and found that there were faster load times when these files existed. It was 6 seconds faster (27.3 sec versus 33.3 sec) compared with/without .sbn and .sbx files.
.cpg are optional plain text files that describes the encoding applied to create the shapefile. If your shapefile doesn’t have a cpg file, then it has the system default encoding.

Mapping Examples

Example 1: Dammed California

Data: California Jurisdictional Dams

Accessed from: https://hub.arcgis.com/datasets/98a09bec89c84681ae1701a2eb62f599_0/data?geometry=-150.074%2C31.096%2C-87.54%2C43.298&page=10

“This dataset is a feature class identifying all dams currently under the jurisdiction of the Division of Safety of Dams (DSOD). The dataset is extracted from DSOD internal records and contains basic information about the dam including the type of construction, basic dimensions such as height, length, and maximum storage capacity; abbreviated owner information to identify the entity legally responsible for the dam; an assessment of the downstream hazard associated with the dam; an assessment of the current condition of the dam; and indication as to whether the dam is operating at a restricted storage level. Several dams span rivers that define county boundaries, so DSOD references the right abutment of the dam to identify the location of the structure and to associate it with a singular administrative subdivision of California.”

Data: California eco-regions (EPA)

Accessed from: https://www.epa.gov/eco-research/ecoregion-download-files-state-region-9

Read in the California ecoregions data (layer “ca_eco”), select only the attribute for eco-region (US_L3NAME), rename that to “Region”, simplify the polygons (for time) using st_simplify, and set the CRS:

Example 1. Dams in CA

#has point location of dams, and other attributes - name, height, ownership, hazars...
ca_eco <- read_sf(dsn = ".", layer = "ca_eco") %>% # Get data! getting simple features information; "." means get from current working directory; then tell it that the layers we want to bring in all start with "ca_eco", not other, really heavy layers
  dplyr::select(US_L3NAME) %>% # Only select column with eco-regions
  rename(Region = US_L3NAME) %>% # Rename that column to "Region"
  st_simplify(dTolerance = 10) %>% # Simplify polygons (for time) - the higher tolerance, the higher resolution
  st_transform(crs = 4326) # Change CRS to 4326

st_crs(ca_eco) #Check projection using st_crs

## Coordinate Reference System:
##   EPSG: 4326 
##   proj4string: "+proj=longlat +datum=WGS84 +no_defs"

Read in the California Counties shapefile data, and set CRS:

ca_counties <- read_sf(dsn = ".", layer = "california_county_shape_file") # Read data
st_crs(ca_counties) # we get N/A - we don't have a CRS set up yet. So make one:

## Coordinate Reference System: NA

st_crs(ca_counties) = 4326 # Set CRS
st_crs(ca_counties) #now we see it has CRS data

## Coordinate Reference System:
##   EPSG: 4326 
##   proj4string: "+proj=longlat +datum=WGS84 +no_defs"

Dams in CA

ca_dams <- read_sf(".", layer = "California_Jurisdictional_Dams") %>% 
  rename(Condition = Condition_) #there was a column with an underscore at the end- I changed the column name/header

ca_dams$Condition <- fct_relevel(ca_dams$Condition, "Fair", "Satisfactory", "Unsatisfactory", "Poor") #rearrange by condition
#***default is alphabetical***, so tell it to give it a different order
levels(ca_dams$Condition)

## [1] "Fair"           "Satisfactory"   "Unsatisfactory" "Poor"

2.0 Maps!

plot(ca_eco)

plot(ca_counties) #autoplots by all indicators

3.0 Make a map with ggplot!

#make a 9-level color brower pallete, and set to 13-levels
#by number of regions
color_count <- 13
#making my own palette based on colorbrewer pallete
my_colors <- colorRampPalette(brewer.pal(10, "Set2"))(color_count) #instead of "10", you can just say "3" and it'll pick the FIRST 3 colors in the  palette and redistribute them by shades

## Warning in brewer.pal(10, "Set2"): n too large, allowed maximum for palette Set2 is 8
## Returning the palette you asked for with that many colors

ggplot(ca_eco) + #baselayer of ecoregion
  geom_sf(aes(fill = Region),
          color = "NA", #color of lines - remove
          show.legend =  FALSE) +  #create our map- spatial information
  scale_fill_manual(values = my_colors) +
  geom_sf(data = ca_counties,
          fill = "NA", 
          color = "gray30", 
          size = 0.1) +  #adding another polygon to the map
  geom_point(data = ca_dams,
             aes(x = Longitude, y = Latitude),
             size = 1, color = "gray10", alpha = 0.5) +
  theme_minimal()

  coord_sf(datum = NA) #this makes the x and y axes with the coordinate disappear

## <ggproto object: Class CoordSf, CoordCartesian, Coord, gg>
##     aspect: function
##     backtransform_range: function
##     clip: on
##     crs: NULL
##     datum: NA
##     default: FALSE
##     distance: function
##     expand: TRUE
##     fixup_graticule_labels: function
##     is_free: function
##     is_linear: function
##     label_axes: list
##     label_graticule: 
##     labels: function
##     limits: list
##     modify_scales: function
##     ndiscr: 100
##     range: function
##     render_axis_h: function
##     render_axis_v: function
##     render_bg: function
##     render_fg: function
##     setup_data: function
##     setup_layout: function
##     setup_panel_params: function
##     setup_params: function
##     transform: function
##     super:  <ggproto object: Class CoordSf, CoordCartesian, Coord, gg>

Example 2. Dams in the Sierra Nevada eco-region

sn <- ca_eco %>% 
  filter(Region == "Sierra Nevada") %>%  #all ca_eco data, only pulling out geometry of SN
  st_join(ca_dams)

## although coordinates are longitude/latitude, st_intersects assumes that they are planar

ggplot(sn) +
  geom_sf(data = ca_counties, fill = "wheat3", color = "NA") + #created a map of CA with all county polygons
  geom_sf(fill = "lemonchiffon4", color = "NA") + #line color = NA
  geom_point(aes(x = Longitude, y = Latitude), size = 0.5, color = "red4")

Example 3. Find intersections between polygons:

sb <- ca_counties %>% 
  filter(NAME == "Santa Barbara")
#only SB county polygons (including islands)

eco_clip <- st_intersection(ca_eco, sb) #clip ca_eco polygon by the bounds of SB polygons

## although coordinates are longitude/latitude, st_intersection assumes that they are planar

## Warning: attribute variables are assumed to be spatially constant
## throughout all geometries

plot(eco_clip)

ggplot(eco_clip) +
  geom_sf(data = ca_counties, 
          fill = "grey90", 
          color = "gray80", 
          size = 0.2) +
  geom_sf(aes(fill = Region), color = "NA") + #eco_clip data is the go to, no need to repeat
  scale_fill_manual(values = c("darkolivegreen2", "darkolivegreen", "gold")) +
  coord_sf(xlim = c(-121, -119), ylim = c(33.5, 35.5)) + #set bounding box size by lat, lon
  geom_point(aes(x = -119.6982, y = 34.4208), size = 2) + #manually add point to downtown SB
  geom_text(x = -119.6982, y = 34.35, label = "Santa Barbara") +
  theme_minimal() +
  theme(legend.position = c(0.5,0.15)) +
  labs(x = "", y = "", title = "Santa Barbara County Eco-Regions")

Example 4. Intro to interactive plots with tmap

map_sb_eco <- tm_shape(eco_clip) +
  tm_fill("Region", palette = "RdPu", alpha = 0.5) +
  tm_shape(ca_counties) +
  tm_borders()
  #tm_fill() to fill polygons with color

tmap_mode("view") #sets tmap to interactive viewing

## tmap mode set to interactive viewing

map_sb_eco

#WHAAAAAAAAT!!!!


#leaflet has all kinds of basemaps
leaflet::providers

## $OpenStreetMap
## [1] "OpenStreetMap"
## 
## $OpenStreetMap.Mapnik
## [1] "OpenStreetMap.Mapnik"
## 
## $OpenStreetMap.BlackAndWhite
## [1] "OpenStreetMap.BlackAndWhite"
## 
## $OpenStreetMap.DE
## [1] "OpenStreetMap.DE"
## 
## $OpenStreetMap.CH
## [1] "OpenStreetMap.CH"
## 
## $OpenStreetMap.France
## [1] "OpenStreetMap.France"
## 
## $OpenStreetMap.HOT
## [1] "OpenStreetMap.HOT"
## 
## $OpenStreetMap.BZH
## [1] "OpenStreetMap.BZH"
## 
## $OpenInfraMap
## [1] "OpenInfraMap"
## 
## $OpenInfraMap.Power
## [1] "OpenInfraMap.Power"
## 
## $OpenInfraMap.Telecom
## [1] "OpenInfraMap.Telecom"
## 
## $OpenInfraMap.Petroleum
## [1] "OpenInfraMap.Petroleum"
## 
## $OpenInfraMap.Water
## [1] "OpenInfraMap.Water"
## 
## $OpenSeaMap
## [1] "OpenSeaMap"
## 
## $OpenPtMap
## [1] "OpenPtMap"
## 
## $OpenTopoMap
## [1] "OpenTopoMap"
## 
## $OpenRailwayMap
## [1] "OpenRailwayMap"
## 
## $OpenFireMap
## [1] "OpenFireMap"
## 
## $SafeCast
## [1] "SafeCast"
## 
## $Thunderforest
## [1] "Thunderforest"
## 
## $Thunderforest.OpenCycleMap
## [1] "Thunderforest.OpenCycleMap"
## 
## $Thunderforest.Transport
## [1] "Thunderforest.Transport"
## 
## $Thunderforest.TransportDark
## [1] "Thunderforest.TransportDark"
## 
## $Thunderforest.SpinalMap
## [1] "Thunderforest.SpinalMap"
## 
## $Thunderforest.Landscape
## [1] "Thunderforest.Landscape"
## 
## $Thunderforest.Outdoors
## [1] "Thunderforest.Outdoors"
## 
## $Thunderforest.Pioneer
## [1] "Thunderforest.Pioneer"
## 
## $OpenMapSurfer
## [1] "OpenMapSurfer"
## 
## $OpenMapSurfer.Roads
## [1] "OpenMapSurfer.Roads"
## 
## $OpenMapSurfer.AdminBounds
## [1] "OpenMapSurfer.AdminBounds"
## 
## $OpenMapSurfer.Grayscale
## [1] "OpenMapSurfer.Grayscale"
## 
## $Hydda
## [1] "Hydda"
## 
## $Hydda.Full
## [1] "Hydda.Full"
## 
## $Hydda.Base
## [1] "Hydda.Base"
## 
## $Hydda.RoadsAndLabels
## [1] "Hydda.RoadsAndLabels"
## 
## $MapBox
## [1] "MapBox"
## 
## $Stamen
## [1] "Stamen"
## 
## $Stamen.Toner
## [1] "Stamen.Toner"
## 
## $Stamen.TonerBackground
## [1] "Stamen.TonerBackground"
## 
## $Stamen.TonerHybrid
## [1] "Stamen.TonerHybrid"
## 
## $Stamen.TonerLines
## [1] "Stamen.TonerLines"
## 
## $Stamen.TonerLabels
## [1] "Stamen.TonerLabels"
## 
## $Stamen.TonerLite
## [1] "Stamen.TonerLite"
## 
## $Stamen.Watercolor
## [1] "Stamen.Watercolor"
## 
## $Stamen.Terrain
## [1] "Stamen.Terrain"
## 
## $Stamen.TerrainBackground
## [1] "Stamen.TerrainBackground"
## 
## $Stamen.TopOSMRelief
## [1] "Stamen.TopOSMRelief"
## 
## $Stamen.TopOSMFeatures
## [1] "Stamen.TopOSMFeatures"
## 
## $Esri
## [1] "Esri"
## 
## $Esri.WorldStreetMap
## [1] "Esri.WorldStreetMap"
## 
## $Esri.DeLorme
## [1] "Esri.DeLorme"
## 
## $Esri.WorldTopoMap
## [1] "Esri.WorldTopoMap"
## 
## $Esri.WorldImagery
## [1] "Esri.WorldImagery"
## 
## $Esri.WorldTerrain
## [1] "Esri.WorldTerrain"
## 
## $Esri.WorldShadedRelief
## [1] "Esri.WorldShadedRelief"
## 
## $Esri.WorldPhysical
## [1] "Esri.WorldPhysical"
## 
## $Esri.OceanBasemap
## [1] "Esri.OceanBasemap"
## 
## $Esri.NatGeoWorldMap
## [1] "Esri.NatGeoWorldMap"
## 
## $Esri.WorldGrayCanvas
## [1] "Esri.WorldGrayCanvas"
## 
## $OpenWeatherMap
## [1] "OpenWeatherMap"
## 
## $OpenWeatherMap.Clouds
## [1] "OpenWeatherMap.Clouds"
## 
## $OpenWeatherMap.CloudsClassic
## [1] "OpenWeatherMap.CloudsClassic"
## 
## $OpenWeatherMap.Precipitation
## [1] "OpenWeatherMap.Precipitation"
## 
## $OpenWeatherMap.PrecipitationClassic
## [1] "OpenWeatherMap.PrecipitationClassic"
## 
## $OpenWeatherMap.Rain
## [1] "OpenWeatherMap.Rain"
## 
## $OpenWeatherMap.RainClassic
## [1] "OpenWeatherMap.RainClassic"
## 
## $OpenWeatherMap.Pressure
## [1] "OpenWeatherMap.Pressure"
## 
## $OpenWeatherMap.PressureContour
## [1] "OpenWeatherMap.PressureContour"
## 
## $OpenWeatherMap.Wind
## [1] "OpenWeatherMap.Wind"
## 
## $OpenWeatherMap.Temperature
## [1] "OpenWeatherMap.Temperature"
## 
## $OpenWeatherMap.Snow
## [1] "OpenWeatherMap.Snow"
## 
## $HERE
## [1] "HERE"
## 
## $HERE.normalDay
## [1] "HERE.normalDay"
## 
## $HERE.normalDayCustom
## [1] "HERE.normalDayCustom"
## 
## $HERE.normalDayGrey
## [1] "HERE.normalDayGrey"
## 
## $HERE.normalDayMobile
## [1] "HERE.normalDayMobile"
## 
## $HERE.normalDayGreyMobile
## [1] "HERE.normalDayGreyMobile"
## 
## $HERE.normalDayTransit
## [1] "HERE.normalDayTransit"
## 
## $HERE.normalDayTransitMobile
## [1] "HERE.normalDayTransitMobile"
## 
## $HERE.normalNight
## [1] "HERE.normalNight"
## 
## $HERE.normalNightMobile
## [1] "HERE.normalNightMobile"
## 
## $HERE.normalNightGrey
## [1] "HERE.normalNightGrey"
## 
## $HERE.normalNightGreyMobile
## [1] "HERE.normalNightGreyMobile"
## 
## $HERE.basicMap
## [1] "HERE.basicMap"
## 
## $HERE.mapLabels
## [1] "HERE.mapLabels"
## 
## $HERE.trafficFlow
## [1] "HERE.trafficFlow"
## 
## $HERE.carnavDayGrey
## [1] "HERE.carnavDayGrey"
## 
## $HERE.hybridDay
## [1] "HERE.hybridDay"
## 
## $HERE.hybridDayMobile
## [1] "HERE.hybridDayMobile"
## 
## $HERE.pedestrianDay
## [1] "HERE.pedestrianDay"
## 
## $HERE.pedestrianNight
## [1] "HERE.pedestrianNight"
## 
## $HERE.satelliteDay
## [1] "HERE.satelliteDay"
## 
## $HERE.terrainDay
## [1] "HERE.terrainDay"
## 
## $HERE.terrainDayMobile
## [1] "HERE.terrainDayMobile"
## 
## $FreeMapSK
## [1] "FreeMapSK"
## 
## $MtbMap
## [1] "MtbMap"
## 
## $CartoDB
## [1] "CartoDB"
## 
## $CartoDB.Positron
## [1] "CartoDB.Positron"
## 
## $CartoDB.PositronNoLabels
## [1] "CartoDB.PositronNoLabels"
## 
## $CartoDB.PositronOnlyLabels
## [1] "CartoDB.PositronOnlyLabels"
## 
## $CartoDB.DarkMatter
## [1] "CartoDB.DarkMatter"
## 
## $CartoDB.DarkMatterNoLabels
## [1] "CartoDB.DarkMatterNoLabels"
## 
## $CartoDB.DarkMatterOnlyLabels
## [1] "CartoDB.DarkMatterOnlyLabels"
## 
## $HikeBike
## [1] "HikeBike"
## 
## $HikeBike.HikeBike
## [1] "HikeBike.HikeBike"
## 
## $HikeBike.HillShading
## [1] "HikeBike.HillShading"
## 
## $BasemapAT
## [1] "BasemapAT"
## 
## $BasemapAT.basemap
## [1] "BasemapAT.basemap"
## 
## $BasemapAT.grau
## [1] "BasemapAT.grau"
## 
## $BasemapAT.overlay
## [1] "BasemapAT.overlay"
## 
## $BasemapAT.highdpi
## [1] "BasemapAT.highdpi"
## 
## $BasemapAT.orthofoto
## [1] "BasemapAT.orthofoto"
## 
## $nlmaps
## [1] "nlmaps"
## 
## $nlmaps.standaard
## [1] "nlmaps.standaard"
## 
## $nlmaps.pastel
## [1] "nlmaps.pastel"
## 
## $nlmaps.grijs
## [1] "nlmaps.grijs"
## 
## $nlmaps.luchtfoto
## [1] "nlmaps.luchtfoto"
## 
## $NASAGIBS
## [1] "NASAGIBS"
## 
## $NASAGIBS.ModisTerraTrueColorCR
## [1] "NASAGIBS.ModisTerraTrueColorCR"
## 
## $NASAGIBS.ModisTerraBands367CR
## [1] "NASAGIBS.ModisTerraBands367CR"
## 
## $NASAGIBS.ViirsEarthAtNight2012
## [1] "NASAGIBS.ViirsEarthAtNight2012"
## 
## $NASAGIBS.ModisTerraLSTDay
## [1] "NASAGIBS.ModisTerraLSTDay"
## 
## $NASAGIBS.ModisTerraSnowCover
## [1] "NASAGIBS.ModisTerraSnowCover"
## 
## $NASAGIBS.ModisTerraAOD
## [1] "NASAGIBS.ModisTerraAOD"
## 
## $NASAGIBS.ModisTerraChlorophyll
## [1] "NASAGIBS.ModisTerraChlorophyll"
## 
## $NLS
## [1] "NLS"
## 
## $JusticeMap
## [1] "JusticeMap"
## 
## $JusticeMap.income
## [1] "JusticeMap.income"
## 
## $JusticeMap.americanIndian
## [1] "JusticeMap.americanIndian"
## 
## $JusticeMap.asian
## [1] "JusticeMap.asian"
## 
## $JusticeMap.black
## [1] "JusticeMap.black"
## 
## $JusticeMap.hispanic
## [1] "JusticeMap.hispanic"
## 
## $JusticeMap.multi
## [1] "JusticeMap.multi"
## 
## $JusticeMap.nonWhite
## [1] "JusticeMap.nonWhite"
## 
## $JusticeMap.white
## [1] "JusticeMap.white"
## 
## $JusticeMap.plurality
## [1] "JusticeMap.plurality"
## 
## $Wikimedia
## [1] "Wikimedia"

#do this to view providers of maps

tm_basemap("CartoDB.DarkMatter") +
  tm_shape(eco_clip) +
  tm_borders(col = "white")

#now no options for basemaps

Example 5. Fault Lines

Fault line data from California Dept. of Conservation: https://maps.conservation.ca.gov/geology/#datalist Separate fault line types syncline/anticline, certain/concealed, direction columns using tidyr::separate().

fault_lines <- read_sf(".", layer = "GMC_str_arc") %>% 
  st_transform(crs = 4326) %>% #use if there is an existing projection; otherwise st_crs() to set it. If the file has a ".prj", it has a projection
  separate(LTYPE, into = c("syn_ant", "certainty", "direction"), sep = ",") #use the comma within the string to separate the words; will populate some with NAs if there were less than 3 arguments within original LTYPE

## Warning: Expected 3 pieces. Missing pieces filled with `NA` in 1091 rows
## [1, 2, 3, 4, 5, 6, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21,
## 22, ...].

plot(fault_lines) #not concerning at all...

ggplot() +
  geom_sf(data = ca_counties, fill = "black", color = "NA") +
  geom_sf(data = fault_lines, aes(color = syn_ant)) +
  theme_void()

# only want to look at faultlines in SB county
sb_faults <- fault_lines %>% 
  st_intersection(sb) #find faultlines only with sb

## although coordinates are longitude/latitude, st_intersection assumes that they are planar

## Warning: attribute variables are assumed to be spatially constant
## throughout all geometries

ggplot() +
  geom_sf(data = sb) +
  geom_sf(data = sb_faults, aes(color = syn_ant))

tm_basemap("CartoDB.DarkMatter") +
  tm_shape(sb) +
  tm_borders(col = "gray50", lwd = 2) +
  tm_shape(sb_faults) +
  tm_lines(col = "syn_ant", palette = c("orange", "purple"), lwd = 2) #you give it the column name in quotations, not within an aes

Example 6. Faceted Maps

ggplot() + 
  geom_sf(data = ca_counties, fill = "black", color = "NA") +
  geom_sf(data = fault_lines, aes(color = syn_ant)) +
  facet_wrap(~syn_ant)

tm_basemap("CartoDB.DarkMatter") +
  tm_shape(sb) +
  tm_borders(col = "gray50", lwd = 2) +
  tm_shape(sb_faults) +
  tm_lines(col = "syn_ant", palette = c("orange","purple"), lwd = 2) +
  tm_facets(by = "syn_ant")

Example 7. Making Spatial Points

ca_sites <- read_csv("cadfw_sensitive_sites.csv")

## Parsed with column specification:
## cols(
##   X = col_double(),
##   Y = col_double(),
##   OBJECTID = col_double(),
##   ACP_NUM = col_double(),
##   SITE_NUM_N = col_character(),
##   SITE_NAME = col_character(),
##   PRI_CODE = col_character(),
##   LATDD = col_double(),
##   LONDD = col_double()
## )

sites_sf <- st_as_sf(ca_sites, coords = c("LONDD", "LATDD"), crs = 4326) #now R knows this is spatial data (otherwise those are just variables and values)
#it took those coordinates and made them sticky geometry!!! 

ggplot() +
  geom_sf(data = ca_counties, fill = "gray40") +
  geom_sf(data = sites_sf, aes(color = PRI_CODE), size = 0.5)   #would have used geom_point for NON-sticky geom data. but now is sticky geometry

Example 8. Choropleth of CA counties by NUMBER of dams in each

#Find counts of dams per county:
intersection <- st_intersection(x = ca_dams, y = ca_counties)

## although coordinates are longitude/latitude, st_intersection assumes that they are planar

## Warning: attribute variables are assumed to be spatially constant
## throughout all geometries

dams_per_county <- intersection %>% 
  group_by(NAME) %>% 
  tally() #will count per name of county

# Check it out: 
# View(dams_per_county)

# Then merge to the ca_counties data: 
ca_tot <- ca_counties %>% 
  st_join(dams_per_county) %>% 
  dplyr::select(NAME.x, n) %>%
  rename(name = NAME.x)

## although coordinates are longitude/latitude, st_intersects assumes that they are planar

# Reassign NA values to zero:
ca_tot$n[is.na(ca_tot$n)] <- 0 #because if there is county without dams, the county itself wont show up. So by assigning 0, keep the county there

#plot
ggplot() +
  geom_sf(data =ca_tot, aes(fill = n), size = 0.2) +
  theme_minimal() +
  scale_fill_continuous(low = "yellow", high = "red")

ESM 244 lab 6

Ronnie Bailey-Steinitz

2/14/2019